#pragma once
#include "elisa/elisa.h"
#include <iostream>
#include <algorithm>
#include <string>
#include <vector>
#include <gaol/gaol>

#include "Modeler.h"

#ifdef ELISA_PERFORMER
#include "elisa_chrono.h"
#endif

using namespace std;
using namespace elisa;

/* Problem: Chemical Equilibrium System

@Article{MorganSommeseWampler1992,
  author    = {A. P. Morgan and A. J. Sommese and C. C. Wampler},
  title     = {{C}omputing {S}ingular {S}olutions to {P}olynomial {S}ystems},
  year      = {1992},
  journal   = {Advances in Applied Mathematics},
  volume    = {13},
  pages     = {305--327}
}
*/
class Chemistry1 : public Modeler
{
private:
	Model _model;

	RealVar z1;
	RealVar z2;
	RealVar z3;
	RealConstraint c1;
	RealConstraint c2;
	RealConstraint c3;

public:
	Chemistry1()
		: z1(Interval(-1000,1000), "z1"),
		  z2(Interval(-1000,1000), "z2"),
		  z3(Interval(-1000,1000), "z3"),
		  c1(14*pow(z1,2) + 6*z1*z2 + 5*z1 - 72*pow(z2,2) - 18*z2 == 850*z3 - 2.0e-9),
		  c2(0.5*z1*pow(z2,2) + 0.01*z1*z2 + 0.13*pow(z2,2) + 0.04*z2 == 4.0e4),
		  c3(0.03*z1*z3 + 0.04*z3 == 850),
		  _model()
	{
	}

	elisa::Model & model()
	{
		_model = c1 && c2 && c3;
		return _model;
	}

	void print()
	{
		cout << '\t' << z1.getName() << " = " << z1.domain() << endl;
		cout << '\t' << z2.getName() << " = " << z2.domain() << endl;
		cout << '\t' << z3.getName() << " = " << z3.domain() << endl;
	}

	RealVar & varX()
	{
		return z1;
	}
	RealVar & varY()
	{
		return z2;
	}
};
